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We present an extensive study of Mott insulator (MI) and superfiuid (SF) shells in Bose-Hubbard 
(BH) models for bosons in optical lattices with harmonic traps. For this we develop an inhomo- 
geneous mean-field theory. Our results for the BH model with one type of spinless bosons agrees 
quantitatively with quantum Monte Carlo (QMC) simulations. Our approach is numerically less 
intensive than such simulations, so we are able to perform calculation on experimentally realistic, 
large 3D systems, explore a wide range of parameter values, and make direct contact with a variety 
of experimental measurements. We also generalize our inhomogeneous mean-field theory to study 
BH models with harmonic traps and (a) two species of bosons or (b) spin-1 bosons. With two 
species of bosons we obtain rich phase diagrams with a variety of SF and MI phases and associated 
shells, when we include a quadratic confining potential. For the spin-1 BH model we show, in a 
representative case, that the system can display alternating shells of polar SF and MI phases; and 
we make interesting predictions for experiments in such systems. 

PACS numbers: 05.30Jp, 67.40Db, 73.43Nq 



I. INTRODUCTION 



High-precision experiments on cold atoms, such as 
spin-polarized 87 Rb, in traps have provided powerful 
methods for the study of quantum phase transitions 1 , 
e.g., the transition from a superfiuid (SF) to a bosonic 
Mott-insulator (MI) in an optical lattice^. This transi- 
tion was predicted by mean-field studies^ and obtained 
in Monte-Carlo simulations 6 of the Bose-Hubbard model 
before it was seen in experiments^—. Recent experi- 
ments 7 -^ have investigated a heteronuclear degenerate 
mixture of two bosonic species, e.g., 87 Rb and 41 K, in 
a three-dimensional optical lattice; such mixtures have 
also been studied theoretically^— and by Monte Carlo 
simulations 14 . Systems of alkali atoms with nuclear spin 
7 = 3/2 have hyperfine spin F = 1; examples include 
23 Na, 39 K, and 87 Rb; these spins are frozen in magnetic 
traps, so these atoms are treated as spinless bosons; how- 
ever, in purely optical traps, such spins can form spinor 
condensates^—. Thus, we consider the following three 
types of Bose-Hubbard (BH) models: (1) a BH model for 
spinless interacting bosons of one type; (2) a generaliza- 
tion of the spinless BH model with two types of bosons; 
and (3) a spin-1 generalization of the spinless BH model 
with bosons of one type. We study these models by de- 
veloping extensions of an inhomogeneous mean-field the- 



ory^, which has been used for the Bose-glass phase in 
the disordered BH model. 

In addition to the optical-lattice potential, a confin- 
ing potential, which is typically quadratic, is present in 
all experiments. This inhomogeneous potential leads to 
inhomogeneities in the phases that are obtained: simu- 
lation o 20 ' 21 of the Bose-Hubbard model with this confin- 
ing potential and experiment s 22 ' 23 on interacting bosons 
in optical lattices with a confining potential have both 
seen alternating shells of SF and MI regions in the single- 
species, spinless case. We explore such shells via the in- 
homogeneous mean-field theory, first for single-species, 
spinless bosons and then for the two-species and spin-1 
generalizations mentioned above. 

Mean-field theories for the Bose-Hubbard model were 
first developed for the homogeneous case 4 ^; these the- 
ories were then extended to the inhomogeneous case 1 ^ 
to develop an understanding of the Bose-glass phase in 
the disordered Bose-Hubbard model. We show that BH 
models with confining potential can be treated, at the 
level of mean-field theory, as was done in the Bose-glass 
case 1 ^; in particular, we provide a natural framework for 
understanding alternating SF and MI shells, which are 
seen in simulation s] 20 ' 21 and experiment er 22 ' 23 on inter- 
acting bosons, trapped in a confining potential, and in 
an optical lattice. Though other groups 2 ^— have stud- 
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ied such shell structure theoretically, they have not ob- 
tained the quantitative agreement with quantum Monte 
Carlo (QMC) simulations^ that we obtain, except in one 
dimension^. Furthermore, our theory yields results in 
good agreement with a variety of experiments; and it 
can be generalized easily to (a) two species of interacting 
bosons and (b) the spin-5 1 case, as we show explicitly for 
S = 1; in both these cases we provide interesting pre- 
dictions that will, we hope, stimulate new experiments. 
Our inhomogeneous mean-field calculations can be car- 
ried out with experimentally realistic parameters, so we 
can make direct comparison with experiments. In par- 
ticular, we obtain in-trap density distributions of alter- 
nating SF and MI shells; these show plateaux in certain 
regions, which can be understood on the basis of sim- 
ple geometrical arguments. Furthermore, we obtain the 
radii of SF and MI shells from in-trap density distribu- 
tions and demonstrate how the phase diagram of the ho- 
mogeneous Bose-Hubbard model can be obtained from 
these radii. We also obtain results that are of direct rel- 
evance to recent atomic-clock-shift experiments^. With 
two species of bosons we obtain phase diagrams in the 
homogeneous case over a far wider range of parameters 
than has been reported hitherto. We find rich phase di- 
agram with phases that include ones in which (a) both 
types of bosons are in SF states, (b) both types of bosons 
are in MI phases with different or the same densities, and 
(c) one type of boson is in an SF phase whereas the other 
type is in an MI phase. We show that each of these phases 
appear in shells when we include a quadratic confining 
potential; and we also obtain in-trap density distribu- 
tions that shows plateaux as in the single-species case. 
In the case of the spin-1 Bose-Hubbard model we show, 
in a representative case, that the system can display al- 
ternating shells of polar SF— and MI phases; the latter 
have integral values for the boson density. Our inhomo- 
geneous mean-field theory leads to interesting predictions 
for atomic-clock-shift experiments in systems with spin-1 
bosons in an optical lattice with a confining potential. 



The remaining part of this paper is organized as fol- 
lows. In Sec. 2 we describe the models we use and how we 
develop an inhomogeneous mean-field theory for them. 
Section 3 is devoted to our results; subsection 3 A con- 
tains the results of our inhomogeneous mean-field the- 
ory for the single-species, spinless Bose-Hubbard model; 
subsection 3B is devoted to the results, for both homo- 
geneous and inhomogeneous cases, for the spinless BH 
model with two types of bosons; subsection 3C is de- 
voted to our results for the single-species BH model for 
spin-1 bosons. Section 4 contains our conclusions, a com- 
parison of our work with earlier studies in this area, and 
the experimental implications of our results. 



II. MODELS AND INHOMOGENEOUS 
MEAN-FIELD THEORY 

We begin by defining the three Bose-Hubbard models 
that we study. We then develop inhomogeneous mean- 
field theories that are well suited for studying the spatial 
organization of phases in these models with confining po- 
tentials. 



A. Models 

The simplest Bose-Hubbard model describes a single 
species of spinless bosons in an optical lattice by the fol- 
lowing Hamiltonian: 



Tt = ~\ E (°<V 



-h. 



Ui (TU— 1) 



<i,j> 



(1) 

here spinless bosons hop between the z nearest-neighbor 
pairs of sites < i,j > with amplitude t, <zj, aj, and 
rii = a\a>i are, respectively, boson creation, annihilation, 
and number operators at the sites i of a d-dimensional 
hypercubic lattice (we study d = 2 and 3), U the onsite 
Hubbard repulsion, fa = (j,— VrRi, fJ> the uniform chem- 
ical potential that controls the total number of bosons, 
Vt the strength of the harmonic confining potential, and 

&i = En=i^n(')i where X n (i), 1 < n < d, arc the 
Cartesian coordinates of the site i (in d = 3 X\ = X, 
X2 = Y, and X3 = Z) ; the origin is chosen to be at the 
center of the lattice. In terms of experimental parame- 
W 



tersi = -^-^fe 2 ^^ 7 , where E r is the recoil energy, 
Vb the strength of the lattice potential, a s (= 5.45nm 
for 87 Rb) the s-wave scattering coefficient, a — A/2 the 
optical lattice constant, and A = 825nm the wavelength 
of the laser used to create the optical lattice; typically 
< Vq < 22E r . We set the scale of energies by using 
zt = 1 in the Bose-Hubbard model [IJ for comparisons 
with experimental systems we should scale all energies 
by E r . 

For a mixture with two types of bosons, we use the 
following Bose-Hubbard Hamiltonian: 



n 

z 



-7 E ^ + M - 7 E ( 6 & + M ( 2 ) 

<i,j> <i,j> 
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the first and second term represent, respectively, the hop- 
ping of bosons of types a and b between the nearest- 
neighbor pairs of sites < i,j > with hopping amplitudes 



t a and tb; here a 



and 



t 



a\a,i and b\, bi, and 



hu — b\bi are, respectively, boson creation, annihilation, 
and number operators at the sites i of a d-dimensional 
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hypercubic lattice (we study d = 2 and 3) for the two 
bosonic species. For simplicity we restrict ourselves to 
the case t a = tb = t, and, to set the scale of energy, we use 
zt = 1. The third and fourth terms account for the onsite 
interactions of bosons of a given type, with energies U a 
and Ub, respectively, whereas the fifth term, with energy 
U a b, arises because of the onsite interactions between 
bosons of types a and b. We have two chemical potential 
terms: p ai = p a -V Ta R^ and pu = Pb~V T bRh where p a 
and \ib the uniform chemical potentials that control the 
total number of bosons, of species a and b, respectively, 
Vra and Vrb are the strengths of their harmonic confining 
potentials (we restrict ourselves to the case Vr a — Vrb), 
and Rf = Yl n =i -^n(*)> where X n (i), 1 < n < d, are the 
Cartesian coordinates of the site i; the origin is at the 
center of the lattice. 

The spin-1 Bose-Hubbard Hamiltonian 1 - that we con- 
sider is 



H 

zt 



<i,j>,cr 



i 

Ifii) -y^/ijnj. (3) 



2 zt ^ { 1 



Here spin-1 bosons can occupy the sites i of a 
d— dimensional, hypercubic lattice and hop between the 
z nearest-neighbor pairs of sites < i,j > with ampli- 
tude t, a is the spin index that can be 1,0,-1, a\ a 
and di^a are, respectively, site- and spin-dependent bo- 
son creation and annihilation operators, and the num- 
ber operator hi a = a[ a ai^', the total number operator 
at site i is hi = ^2 a hi, a , and Fi = J2a,a> a\ a F^ a , a^* 

with Fc^i standard spin-1 matrices. The model ([3]) in- 
cludes, in addition to the onsite repulsion Uq, an energy 
U2 for nonzero spin configurations on a site. Such a spin- 
dependent term arises from the difference between the 
scattering lengths for S = and S = 2 channels 31 . The 
inhomogeneous chemical potential pi is related to the 
uniform chemical potential p and the quadratic, confin- 
ing potential as in the spinless case[TJ We set the scale 
of energies by choosing zt = 1. 



B. Inhomogeneous mean-field theory 

The mean-field theory we use has been very successful 
in obtaining the phase diagrams for models ((T|) and 
with Vt — 0, i.e., in the absence of the harmonic confin- 
ing potentia&i£. The inhomogeneous generalization of 
this theory, developed first for the Bose-glass phased and 
spinless bosons, decouples the hopping term to obtain an 
effective one-site problem, neglects quadratic deviations 
from equilibrium values (denoted by angular brackets), 
uses the approximation 



aja. 



(aha* + a] (a 



introduces the superfluid order parameter ipi = (ai) for 
the site i, and thence expresses the Hamiltonian (TTJ) as 
H MF = J2i'R-i IF '■> where the superscript MF denotes 
mean field and the single-site Hamiltonian is 



MF 



Zt 



1 U 

Y7t 



hi(h i -l)-^h i -((f) l al + 4 l *a i )+-ip*(l) i . (5) 



zt 



Here (pi = | J2s V'i+i and 5 labels the z nearest neigh- 
bors of the site i. If Vr = 0, the effective onsite chemical 
potential pi = p, for all i, so the local density and su- 
perfluid order parameters are independent of i: pi = p 
and i/ji = ip. If Vt > we first obtain the matrix el- 
ements of Hf IF in the onsite, occupation-number basis 
{|rij)}, truncated in practice by choosing a finite value for 
Umax, the total number of bosons per site, for a given 
initial set of values for {ipi}. [For small values of U we 
must use large values of nmax; for the values of U we 
consider n ma x = 6 suffices.] We then diagonalize this 
matrix, which depends on ipi and ipi+$, to obtain the 
lowest energy and the corresponding wave function, de- 
noted, respectively, by E g (tpi,tpi^s) and 5 , ff ({-0i}); from 
these we obtain the new superfluid order parameters 

ipi = i^gd^i}) I a i I ^'^({V'i})); we use these new values 
of ipi as inputs to reconstruct Hf IF and repeat the diag- 
onalization procedure until we achieve self consistency of 
input and output values to obtain the equilibrium value 
ipl q (henceforth we suppress the superscript eq for nota- 
tional convenience) . [This is equivalent to a minimization 
of the total energy E g ({ipi}) = J2i ^K^i^i+s) with re- 
spect to ipi] if more than one solution is obtained, we pick 
the one that yields the global minimum.] The onsite den- 
sity is obtained from p, = {9 g ({ipi}) | hi | ^({V'i}))- In 
representative cases, we have found that the equilibrium 
value of ipi is real; so, henceforth, we restrict ourselves to 
real values of ip>. 

For the two-species Hamiltonian ([3]) our mean-field 
theory obtains an effective one-site problem by decou- 
pling the two hopping terms as follows (cf., Eq. [4j: 



4 a j - ( a l) a j + a \ ( a j ) - (4 ){aj); 

b\b 3 ~ + -$><&,■>; (6) 



here the superfluid order parameters for the site i for 
bosons of types a and b are ip a i = (ai) and ipbi = (f>i)i 
respectively. The approximation ([5]) can now be used 
to write the Hamiltonian |3) as a sum over single-site, 
mean-field Hamiltonians Hf IF (cf. , Eq. [SJ given below: 



MF 



Zt 



1 Ua / .* -. \ Pai ~ 

77 T^aiy^ai T^az 

2 zt zt 

-{4>ar4 J r(p* ai a l ) + V>a Ai 

1 Ub A t * > nu „ 
+ -—nbi{nu - 1) -nu 

2 zt zt 



(7) 



kib\ + (pbih) + tpbi4>bi + ^-haihb 



4)(a 3 ) 



(4) 



Here (f> ai = \ J2s and cp u = \ J^s ipbi+6, where S 

labels the nearest neighbors of the site i. If Vt = 0, the 
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effective onsite chemicai potentials p, a i = p a and pu = 
fi b , for all i, so p ai = p a , p bi = p b , ip ai = ip a , and ip bl = ip b 
are independent of i. 

If Vt > 0, we first obtain, for a given initial set of val- 
ues for {ipai} and {ipbi}, the matrix elements of Wf' IF 
in the onsite, occupation-number basis {|n a j) ,\nu)}, 
which we truncate in a practical calculation by choos- 
ing a finite value rimax for the total number of bosons 
per site. [The smaller the values of the interaction pa- 
rameters U a ,Ub, and U a b the larger must be the value 
of n max ; for the values of U a , U b , U ab , p a , and p b 
that we consider, nmax = 6 suffices.] We then di- 
agonalize this matrix, which depends on ip a i , ipbi, 
4>a{i+8) , and tpbu+s) to obtain the lowest energy and the 
corresponding wave function, denoted, respectively, by 

Eg(.i>ai,ip a (i+6); ipbi,i>b(i+5)) and i &g({ipai,i>bi}), whence 
we obtain the new superfiuid order parameters ip a i = 
{^ g ({ipai,ipbi}) | a-i | ^ g {{>Pa,i,ipu})) and ipu = 
(^g({ip a i,^bi}) I bi | $g({ip a i,ip b i})); we use these new 
values of ip a i and ipu as inputs to reconstruct Wf IF and 
repeat the diagonalization procedure until we achieve self 
consistency of input and output values to obtain the 
equilibrium value ip^i and iplf; again we suppress the 
superscript eq for notational convenience. [As we have 
mentioned in the single-species case, this self-consistency 
procedure is equivalent to a minimization, with respect 
to ipai and ipu, of the total energy E g ({ip ai ,ip b i}) = 
J2i E g (ipai,ipai+s; ipu,ipbi+s)] we pick the one that yields 
the global minimum.] The onsite densities are obtained 
from pai = {9 g ({ip a i,ipu}) I n ai | ^ g ({ipai, ipbi})) and 
Pbi = (^g({ipai,ipbi}) | flu | ^ g ({ipai, ipbi})}, respectively. 
We follow our discussion of the mean-field theory of the 
BH model (TIJ and restrict ourselves to real values of ip a i 
and ip bi . 

The inhomogeneous mean-field theory for the spin-1 
BH model follows along similar lines. The spin-1 analogs 
of Eqs. S] and [S] are respectively, 



and 

— - — = nArii — 1) H [r, — 2rii) n; 

zt 2 zt V ' 2 zt V 1 ' zt 

- ^2(<t>i,<jat,o- + <Pi* a i,<?) + ^tv&hv ( 9 ) 

cr cr 

Here we use the following superfiuid order parameters: 



4>i,a 



(10) 



and 4>i,a = h E<5 ^(i+<S),<ri where and 5 labels the z near- 
est neighbors of the site i; recall, furthermore, that a 
can assume the values 1, 0, —1, and h ia = a\ a a ii(J , 

™» = E CT ™i,<n and = J2cr,a' a i, a F*,<T' a i,<r' with F<r,<r' 

standard spin-1 matrices. With these order parame- 
ters (cf., Eq. [TU]) we have developed an inhomogeneous 
version of the homogeneous mean-field theory^ for the 
spin-1 BH model with Vt = 0. 



The self-consistency procedure that we use now is sim- 
ilar to, but more complicated than, the one we have used 
for the spinless BH model. If Vt > we first obtain, 
for a given initial set of values for {ipi, a }, the matrix el- 
ements of Hf IF in the onsite, occupation-number basis 
{In^-i, fiifi, rii t i)}, truncated in a practical calculation 
by choosing a finite value for rimax, the total number 
of bosons per site, [For small values of U and U2 we 
must use large values of nmax! for the values we use 
here, rimax = 4 suffices.] We then diagonalize this ma- 
trix, which depends on ip^ a and ip^+s),a, to obtain the 
lowest energy and the corresponding wave function, de- 
noted, respectively, by E l g (ip i>(r , ip^ +S ) j(T ) and ^ g ({ipi,a}); 
from these we obtain the new superfiuid order parame- 
ters ip it<T = I o,> I ^ g ({ipi,a})); we use these 
new values of ipi^ a to reconstruct Hf ,IF and repeat the 
diagonalization procedure until input and output values 
are self consistent; thus we obtain the equilibrium value 
V'iff- We suppress eq as above and recall that this self- 
consistent procedure is equivalent to a minimization of 
the total energy in the spin-1 case 18 . Here too, we fol- 
low our discussion of the mean-field theory of the BH 
model H]) and restrict ourselves to real values of tpi !<7 . 

We have noted in an earlier studyi^ that, at the level 
of our mean-field theory, the superfiuid density in the 
spin-1 case is 



= £ic 9 i 2 ; 



(ii) 



and the magnetic properties of the SF phases follow 
fromikil 



e <?|2 



(12) 



If we substitute the explicit forms of the spin-1 matrices 
we obtain 

{ f) = ^t+h^t+^-f-?*, 



(F) 1 



i W>iVo + V'-iV'o) 2 , (iPl-iP 2 -i) 2 



, (13) 



where x and z are unit vectors in spin space; SF phases 
with (F) = and (F) 2 = 1 are referred to as polar and 
ferromagnetic, respectively. The order-parameter mani- 
folds of these phases can be found in earlier studie o 17 ' 18 . 



III. RESULTS 

Given the formalism we have developed above, we can 
obtain several results for quantities that have been mea- 
sured in quantum Monte Carlo (QMC) simulations or in 
experiments for the spinless case. We cover this in Sub- 
section 3A. Subsection 3B is devoted to the results of our 
inhomogeneous MF theory for the case with two types of 
bosons. Subsection 3C is devoted to the results of our 
inhomogeneous MF theory for the spin-1 case. 
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FIG. 1: (Color online) Plots with a comparison of our MF 
results for the local density m (blue filled circles), local com- 
pressibility K l ° cal = dn/dp eff (blue open circles), and local 
superfiuid density pf = ipf (green filled triangles) of spin- 
less bosons in a two-dimensional parabolic trap with Vt/U = 
0.002, n/U = 0.37 and U/t = 25; we have obtained QMC 
data by digitizing plots in figures in simulation studies^ for 
n; (red filled squares) and K l ° cal (red open squares) . 
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FIG. 2: (Color online) (a) SF (white) and MI regions [p = 2 
(red) and p — 1 (black)] annuli formed in a 2D planar section 
V z through the 3D lattice, at a vertical distance z from the 
center (here z = 0) and (b) the corresponding radial variation 
of density pi (red circles) and superfiuid density pf (black 
squares) for p/E r = I, Vr/E r = 0.0003 and Vo = !5E r ; the 
outermost gray regions in (a) contain no bosons. 



A. Results for the spinless Bose-Hubbard model 

First we compare our mean-field (MF) results with 
those obtained by quantum-Monte-Carlo (QMC) sim- 
ulations in two dimensions 20 . These simulations use 
U/t = 25, fi/U = 0.37, and V T /U = 0.002, and ob- 
tain the local density pi and the local compressibility 



K l f cal — dpi/dpLi. For this set of parameters we calculate 
Pi and thence K l ° cal - ) W e also obtain the local superfiuid 
density pf = vpf (the last formula is valid at the level 
of our MF theory). In Fig. [1] we plot versus Ri our MF 
results for pi, n l ° cal , and pf along with data from QMC 
simulations^ -. For this set of parameters, the central re- 
gion near the origin of the lattice is in the MI phase, 
i.e., the local density pi = 1 and both pf and k° van- 
ish. This central core is enveloped by an SF shell, with 
nonzero values for pf and n l ° cal . As we move radially 
outward from the center, pi decreases monotonically till 
it goes to zero, as do pf and K l ° cal , in the region where 
A*i < 0. The quantitative agreement between our MF 
results and those from QMC is shown in Fig. [TJ there 
is only a slight discrepancy between the MF pi and its 
QMC analog at the MI-SF interface; our result for K l P cal 
also seems to miss, at this interface, the shoulder that ap- 
pears in the QMC n l ° cal perhaps because our MF theory 
overestimates the stability of the SF phase. 

This good agreement between our MF results and those 
of QMC simulations has encouraged us to use our MF 
theory in cases where such simulations pose a significant 
numerical challenge. In particular, we use our theory to 
make direct comparisons with experiments^ 2 - that have 
observed alternating MI and SF shells in 3D optical lat- 
tices by recording in-trap density distributions of bosons 
at different filling fractions. We use a simple-cubic lat- 
tice with 121 3 sites, p/E r = 1, V T /E r = 0.0003, and 
the optical potential Vo/E r in the range 12 — 16 so that 
the number of bosons N ~ 10 6 , which is comparable to 
the number of atoms in the experiments^ 2 - we consider. 
This choice of parameters leads to two well-developed MI 
shells (p = 1 and 2, respectively). The MI and SF shells 
appear as annuli 22 in a 2D planar section V z through 
the 3D lattice, at a vertical distance z from the center 
[see, e.g., Fig. [2ja) for Vo/E r = 15 and z — where the 
core region is in the SF phase]. Figure EJb) shows that, 
as we move radially outward, pi decreases monotonically 
and pf is zero in the two MI regimes (16 < Ri < 34 
and 44 < Ri < 52) in which pi is pinned at 2 and 1, 
respectively. SF and MI shells alternate and the outer- 
most one is always in the SF phase; their positions and 
radii depend on p, which also controls the total number 
N of atoms in the system, as illustrated by the V z —q sec- 
tions in Figs[3Ja) and (b) for V /E r = 15 and p = 0.8 
(iV = 7.1 x 10 5 ) and p = 0.9 (N = 8.9 x 10 5 ), respectively. 
For any 2D planar section V z we can calculate N m (z), the 
number of bosons in the p — m MI annulus, and N^z), 
the remaining number of bosons; the total number of 
bosons in this planar section is N(z) = N m (z) + N^z), 
which does not depend on m. In FigsJ3] (c) and (d) we 
show, for m — 2 and p — 0.8 and 0.9, respectively, plots 
versus z of N m (z) (full red squares), N^(z) (full blue tri- 
angles), and their sum N(z) (full black circles). Figures^] 
(c) and (d) are remarkably similar to the density profiles 
obtained in experiments 2 ^ [cf., their Figs. 3(c) and (d)]. 

The radii of MI shells follow from such in-trap den- 
sity profiles: In Figs. BJa) and (b) we plot N m (z) and 
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FIG. 3: (Color online) SF (white) and MI regions [p = 2 (red) 
and p — l(black)] annuli formed in the 2D V z =o sections 
in FigsOa) and (b) for Vo/E r = 15, V T /E r = 0.0003, 
and (a) p/E r = 0.8 (TV = 7.1 x 10 5 ) and (b) p/E r = 0.9 
(TV = 8.9 x 10 5 ). The corresponding integrated in-trap den- 
sity profiles N2(z) (red squares), N^z) (blue squares) and 
7V2 + N2 (black circles) are shown, respectively, in (c) and 
(d); these figures are qualitatively similar to Figs. 3(c) and 
(d) in recent experimental stud}". 



A^(z) versus z for to = 2 and m = 1, respectively, with 
p/E r = 1. The curves N m (z) show nearly flat plateaux 
for —Rj(m) < z < Ri(m); similar plateaux occur in 
N^(z) for Rj(m) <| z \< R (m) [Figs. Ha), (c) and 
(b), (d) for to = 2 and m = 1, respectively]. Here Rj(m) 
and Ro(m) are the inner and outer radii of the MI shell 
with integer density to. Elementary geometry can be 
used to surmise the existence of these plateaux from the 
MI-SF shell structure^ as we show below. 

N m {z) is to times the total number of sites inside the 
p = 777 MI annulus; this number of sites is well approxi- 
mated by the area A(z,m) of this annulus. Thus, 

N m (z) — mA(z, to) = mir{RQ(z, to) — i?f (z, m)], (14) 

where Ro{z, to) and Ri(z, to) are, respectively, the outer 
and inner radii of the MI annulus with density p = to, in 
the 2D planar section V z . If z < Ri (to) , simple geometry 



yields Rj(z,m) = Rj(m) 



therefore, 



and Rq 



(z, to) 



Rl{m) 



N m {z) = mir[R 2 (m) - Rj(m)}, 



(15) 



whence we conclude that N m (z) is independent of z when 
\z\ < Ri(m); this result yields the plateaux in the in-trap 
density profiles shown in Figs. 2] (a)- (d); if \z\ > Ro(m), 



the 2D planar section has no MI shell with density to, 
thus, N m (z) = 0, which is also apparent in these fig- 
ures. For z < Ri(m), the central parts of the 2D planar 
sections V z show SF shells and Nf n = N(z) — N m (z) de- 
creases as we increase z. For Ri(m) < z < Ro{m), the 
central parts of the 2D planar sections V z show MI shells; 
the number of bosons in such MI shells is N m (z) and it 
is proportional to the area of this central shell, namely, 
TO7r(i?Q(TO) — z 2 ); thus, N m (z) decreases as we increase 
z here; however, N 7 m — N(z) — N m {z) remains indepen- 
dent of z, because of the simple geometrical arguments 
given above; i.e., we have plateaux in N 7 m (z) in the region 
Ri{m) < z < R {m). Finally, for z > R {m), the 2D 
planar section V z has no MI shell with density to, from 
which it follows that N m (z) = 0. 
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FIG. 4: (Color online) (a) and (b) Plots of N m (z)(red squares) 
and (blue circles) versus z for m — 2 and 1, respectively, 
with p/E r = 1; the plots in (c) and (d), which are the same 
as the right halves of (a) and (b), respectively, show how 
we determine the inner and outer radii (Ri(m) and Roim), 
respectively) from the plateaux in N m (z) and N^iz). 

In Fig. Eta) we plot, for to = 1 and 2, Ro(m) and 
Ri(m), which we have determined from plots such as 
those in Figs. 0|c) and (d), versus Vo(E r ); the MI phase 
with p = to lies between the curves i?o( m ) an d Ri(m). 
Figure [5ja) can be used to obtain the phase diagram 
of the homogeneous Bose-Hubbard model as follows: 
Pi = fi — VrRi, so Ro{m) and Rj(m) can be used to ob- 
tain /j," (to) = p — VrRoi 171 ) an d p + (m) — p — Vt-Rj(to), 
which are, respectively, the lower and upper boundaries 
of the Mott lobe with density p — m. The resulting Mott 
lobe (obtained by the conversion Vo(E r ) —> U / zt) is given 
in Fig. EJb) along with its counterpart for the homo- 
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FIG. 5: (Color online) (a) Plots, for p = m = 1 and 2, of 
Roijn) and Ri(m) versus Vb(-Er); the MI phase with p — m 
lies between the curves Ro(m) and Ri(m); red circles and red 
inverted triangles denote i?o(l) and Ro(2), respectively; and 
black squares and black triangles denote Ri(l) and Ri(2), 
respectively; (b) Mott-insulating lobes, in plots of p + (m) = 
p — VrRj(m) and p~ (m) = p — Vr-Ro( m ) f° r nn — 2 and 
m — 1 versus U/(zt) (obtained by the conversion Vo(E r ) — > 
U/zt); we use the same symbols as in (a); and we show, for 
comparison, the boundaries of the MI lobes for m = 1 (blue 
diamonds) and m — 2 (blue triangles) that follow from our 
mean-field theory for the homogeneous Bose-Hubbard model 5 . 
(c) An illustrative plot of the total number of bosons N in the 
system versus the chemical potential p. 



geneous Bose-Hubbard model, which we have obtained 
from the homogeneous mean-field theory 5 -; the agreement 
between these lobes is striking; and it encourages us to 
suggest that the phase diagram of the homogeneous Bose- 
Hubbard model can be obtained from the inner and outer 
radii of the MI shells. Thus, experiments on cold atoms in 
optical lattices with a quadratic confining potential, can 
be used directly to obtain the phase diagram of the homo- 
geneous Bose-Hubbard model from Ro(m) and Rj(m), 
which can be determined for an MI shell with density 
p = m as described above. Note that (a) p~(m) and 
p + (m) are fixed for a given Vo(E r ) and (b) the total 
number of bosons N increases linearly with the chemi- 
cal potential p (see Fig. EJc)). Therefore, the inner and 
outer radii of the MI shell Ro,i(m) — \/{p — p~- + /Vr) 
are proportional to for fixed Vr and Vo; this propor- 
tionality has been reported in the recent experiments 22 
[cf., their Fig. 3]. 

Images of MI shells have been obtained recently from 
atomic-clock-shift experiments 2 ^. By using the density- 
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FIG. 6: (Color online) Representative plots of Nb(p), the 
number of bosons in the system with a given density p, versus 
p near (a) p = 1 and (b) p = 2 for Vo = 12E r (black squares), 
14E r (red diamonds) and 16E r (blue triangles), (c) the radial 
variation of density pi for Vo = 12£V (black squares), 14E r 
(red diamonds) and 16E r (blue triangles). 



dependent transition-frequency shifts, sites with differ- 
ent densities of bosons can be distinguished spectroscop- 
ically; and, therefore, MI shells, with different values of 
the integer density m, are revealed as peaks in the oc- 
cupation number at the corresponding frequencies. This 
experiment gives Nb(p), the number of bosons in the sys- 
tem at a given density p. We use our inhomogeneous 
mean-field theory to obtain Nf,(p) and in Figs. |H£a)-(b) 
we plot Nb{p) (with p close to m = 1 and 2, respec- 
tively) for V /E r = 12, 14 and 16, with V T /E r = 0.0003 
and p/E r = 1. The SF and MI shell structure is evident 
from the radial variation of the local density given in Fig. 
\§lc). For Vo — 12E r , no Mott shells is developed; this is 
reflected in a flat variation of Nb(p) for all p. However, 
if Vq = lAE r , there is a well-formed MI shell p = 1; this 
can be inferred from the peak in Nb(p) at p = 1; and, as 
Vq increases, more Mott shells, with higher, integral val- 
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ues of p, appear. This behavior of Nb(p) is in accordance 
with recent experiments^ [cf., their Fig. 1]. 



B. Results for the Bose-Hubbard model with two 
species of bosons 

We begin with an investigation of representative phase 
diagrams of the Bose-Hubbard model ([3]), with two 
species of bosons, in the homogeneous case, i.e., with 
Vto, = Vrb = Vt = 0. These have been explored to some 
extent in earlier theoretical studies^— and Monte Carlo 
simulations 1 ^, but not over as wide a range of parame- 
ters as we consider here. Next we use the inhomogeneous 
mean-held theory that we have developed above to ex- 
plore order-parameter profiles and a variety of MI and 
SF shells that are obtained when we have a quadratic 
trap potential. We also present Fourier transforms of 
one-dimensional sections of these profiles. 

First we consider the case U a b < U a = Ub and p a = 
fib = /i in which the order parameters and densities for 
both types of bosons show the same dependence on p. 
In the first row of Fig. [Jj we show the phase diagram 
(Fig. [JJ (a)), and plots versus p of the order-parameters 
tp a (red line) and ipb (blue dashed line) and the densities 
p a (green dashed line) and pb (pink full line) for U a b = 
0.5U a , U a = U b , and p a = Pb = P and U a = 9 (Fig. [Jj 
(b)), U a = 11 (Fig. [3 (c)), and U a = 13 (Fig. [7J (d)). (We 
do not divide explicitly by zt because we set zt = 1). 
The phase diagram shows an SF phase in which both 
species are superfluid; the blue Mil lobe denotes a Mott- 
insulating phase in which the density p = 1 is attained 
by having p a = Pb = 1/2; the brown MI tt lMI;,l lobe 
denotes a Mott-insulating phase in which the densities 
Pa = Pb = 1; the pink MI a 2MI b 2 lobe denotes a Mott- 
insulating phase in which the densities p a — pt = 2. Such 
phase diagrams can be obtained from plots like those in 
Figs.[7j(b)-(d). 

In the second row of Fig. 1 we show the phase diagram 
(Fig. 1 (e)), and plots versus p of the order-parameters 
ip a and ipb and the densities p a and pb for U a b — 0.2f7 o , 
U b = 0.9U a , and p a = Pb = P and U a = 9 (Fig. [7J (f)), 
U a = 11 (Fig.[7J(g)), and U a = 13 (Fig.[7J(h)). The phase 
diagram shows an SF phase and brown MI a lMI b l and 
pink MI Q 2MI(,2 lobes; these are like their counterparts in 
Fig. [7J(a). In addition we have the following phases: (i) a 
green sliver MI a l in which bosons of type a are in an MI 
phase with p a — 1 and bosons of type b are superfluid; 
(ii) a green-ochre region MI b 2 in which bosons of type b 
are in an MI phase with pb = 2 and bosons of type a are 
superfluid; and (hi) a dark-green region MI 2 in which 
bosons of type a are in an MI phase with p a = 2 and 
bosons of type b are superfluid. 

In the third row of Fig. [7J we show the phase diagram 
(Fig. [7J (i)), and plots versus p of the order-parameters 
ip a and ipb and the densities p a and pb for U a b = 0.6U a , 
U b = O.W a , and p a = Pb = P and U a = 9 (Fig. [7J (j)), 
U a = 11 (Fig. m (k)), and U a = 13 (Fig. 1 (1)). The 



phase diagram shows the following: an SF phase; blue 
Mil, brown MI a l,MI 6 l, and pink MI a 2MI 6 2 lobes; green 
MI a l and dark-green MI a 2 regions; these are like their 
counterparts in Figs. [7J (a) and (e). In addition we have 
a red MI b 2MI a l lobe in which pb = 2 and p a = 1. 

Next we consider the case U a b < U a , Ub — U a , and 
Pb = Q.75p a . Specifically, in the first row of Fig. [8] we 
show the phase diagram (Fig. 1 (a)), and plots versus p a 
of the order-parameters ip a (red line) and ipb (blue dashed 
line) and the densities p a (green full line) and pb (pink 
dashed line) for U a b = 0.1U a , U a = Ub, and pb = 0.75p a 
and U a = 9 (Fig. [5](b)), U a = 11 (Fig. 1(c)), and U a = 13 
(Fig. [5] (d)). (We do not divide explicitly by zt because 
we set zt = 1). The phase diagram shows an SF phase 
and brown MI a lMI h l, pink MI a 2MI fc 2 and red MI a 2MI 6 l 
lobes; green MI a l, dark-green MI a 2 and a green-ochre 
MI b 2 regions; these are like their counterparts in Figs. [7J 
(a) and (e). In addition we have a light-blue region MI b l 
in which bosons of type b are in an MI phase with pb = 1 
and bosons of type a are superfluid. 

In the second row of Fig. [5] we show the phase diagram 
(Fig. [S] (e)), and plots versus p a of the order-parameters 
ip a and ipb and the densities p a and pb for U a b = 0.3C/ a , 
U b = U a , and p b - 0.75p a and U a = 9 (Fig. E (f)), 
U a = U (Fig. M(g)), and U a = 13 (Fig. 1(h)). This 
phase diagram shows the following: an SF phase; brown 
MI Q l,MI fc l and red MI a 2MI b l lobes; green MI a l, dark- 
green MI a 2 and light- blue MIf,l regions; these are like 
their counterparts in Figs. [5] (a). In addition we have a 
dark-gray SF Q phase in which bosons of type a are in an 
SF phase and the bosons of type b have vanished. 

In the third row of Fig. [5] we show the phase diagram 
(Fig. [8] (i)), and plots versus p a of the order-parameters 
tp a and ipb and the densities p a and pb for U a b = 0.7C/ a , 
U b = U a , and p b = 0.75p a and U a = 9 (Fig. 1 (j)), 
U a = 11 (Fig. M (k)), and U a = 13 (Fig. E (!))• The 
phase diagram shows the following: an SF a phase; MI a l 
and MI a 2 regions; these are like their counterparts in 
Figs. E (a) and (e) in which the bosons density for type 
b has vanished. 

We now consider the case U a b < U a , Ub = 0.9U a , and 
Pb = 0-9p a - Specifically, in the first row of Fig.[9]we show 
the phase diagram (Fig. [5] (a)), and plots versus p a of the 
order-parameters %p a (red line) and ipb (blue dashed line) 
and the densities p a (green full line) and pb (pink dashed 
line) for U a b = 0.2[/ a , U a = 0.9C4, and pb = Q-9p a and 
U a = 9 (Fig. 1(b)), U a = 11 (Fig. El ( c)), and U a = 13 
(Fig. 1(d)). The phase diagram shows an SF phase and 
brown MI a lMI fc l and pink MI a 2MI b 2 lobes; green MI Q 1 
and dark-green MI a 2; these are like their counterparts in 
Figs. land Figs. 1(a) and (e). 

In the second row of Fig. 1 we show the phase diagram 
(Fig. 1(e)), and plots versus p a of the order-parameters 
ip a and ipb and the densities p a and pb for U a b = . 5 C/ a , 
U b = 0.9U a , and p b = Q.9p a and U a = 9 (Fig. 1 (f)), 
U a = 11 (Fig. 1 (g)), and U a = 13 (Fig. 1 (h)). The 
phase diagram shows the following: an SF phase and SF a ; 
brown MI a l,MI h l, red MI a 2MI h l and pink MI a 2MI b 2 
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FIG. 7: (Color online) Phase diagram (a), and plots versus p of the order-parameters ip a (red line) and ipb (blue dashed line) 
and the densities p a (green dashed line) and pb (pink full line) for U a b = 0.5U a , U a — Ub, and p a — pb = P and U a = 9 (b), 
U a = 11 (c), and U a = 13 (d). (We do not divide explicitly by zt because we set zt — 1). The phase diagram shows an SF phase 
in which both species are superfluid; the blue Mil lobe denotes a Mott-insulating phase in which the density p = 1 is attained 
by having p a — pb — 1/2; the brown MI a lMI(,l lobe denotes a Mott-insulating phase in which the densities p a — Pb = 1; the 
pink MI a 2MIj,2 lobe denotes a Mott-insulating phase in which the densities p a = pb — 2. In the second row we show the phase 
diagram (e), and plots versus p of the order-parameters ip a and ipb and the densities p a and pb for U a b = 0.2U a , Ub = 0.9J7 O) 
and p a — Pb = p and U a = 9 (f), U a = 11 (g), and U a = 13 (h). The phase diagram shows an SF phase and brown MI a lMIbl 
and pink MI a 2MIi,2 lobes; these are like their counterparts in (a). In addition we have the following phases: (i) a green sliver 
MI a l in which bosons of type a are in an MI phase with p a = 1 and bosons of type b are superfluid; (ii) a green-ochre region 
MI(,2 in which bosons of type 6 are in an MI phase with pb — 2 and bosons of type a are superfluid; and (iii) a dark-green 
region MI a 2 in which bosons of type a are in an MI phase with p a = 2 and bosons of type 6 are superfluid. In the third row we 
show the phase diagram (i), and plots versus p of the order-parameters ip a and ipb and the densities p a and pb for U a b = 0-6U a , 
Ub = 0.9U a , and p a = Pb = P and U a — 9 (j), U a = 11 (k), and U a = 13 (1). The phase diagram shows the following: an SF 
phase; blue Mil, brown MI a l,MI(,l, and pink MI a 2MIi,2 lobes; green MI a l and dark-green MI a 2 regions; these are like their 
counterparts in (a) and (e). In addition we have a red MIt,2MI a l lobe in which pb = 2 and p a = 1. 



lobes; green MI a l, dark-green MI a 2 and light-blue MI^l 
regions; these are like their counterparts in Figs. [S] (a) 
and (e). 

We now consider the effect of a parabolic potential 
and use the inhomogeneous mean-field theory, developed 
in the previous Section, to obtain alternating spherical 
shells of the variety of MI and SF phases, shown in 
the phase diagrams in Figs. [71 [51 and [HI for the two- 
species BH model ([3]). We do this by obtaining the order- 
parameter profiles {ipai, Paii ipbii Pbi} and also by obtain- 
ing in-trap density distributions of bosons at represen- 
tative values of U a b, Ua, Ub, and p a = Pb- In particu- 
lar, we use a 3D simple-cubic lattice with 128 3 sites and 
Vr/{zt) — 0.008; and we study the following represen- 
tative case: p a /(zt) = pbj{zt) = 30 , U a b = 0.6C/ a , 
Ub = 0.9U a , when U a /(zt) = 13. With these parameters 



the total number of bosons Nt — 10 6 , which is compa- 
rable to experimental values. Furthermore, this choice 
of parameters leads not only to SF shells but also two 
well-developed MI shells (Mil and MI2) . 

We show representative plots of the densities p a (green 
dashed line) and pb (pink line) versus the position X 
along the line Y = Z = are shown in Figs. [TU] for 
Pa = Pb = 30, Vr/(zt) = 0.008, and the following six pa- 
rameter sets, respectively: (a) U a /(zt) = 8, Ub — Q-9f7 a , 
and U ab = 0.W a (b) U a /(zt) = 10, U b = 0.9U a , and 
U ab = 0.6U a , (c) U a /(zt) = 11, U b = 0.9U a , and U ab = 
0.6C/ a , (d) U a /(zt) = 12, U b = 0.9U a , and U ab = 0.6£/ a , 
(c) U a /(zt) = 13, U b = 0.9U a , and U ab = 0.6U a , and (f) 
U a /(zt) = 14, U b = 0.9U a , and U ab = 0.6U a . It is also 
useful to obtain a complementary, Fourier-representation 
picture of the profiles in Figs. [TO] (a)- (f) because it might 
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FIG. 8: (Color online) Phase diagram (a), and plots versus fi a of the order-parameters ifj a (red line) and tpt, (blue dashed line) 
and the densities p a (green full line) and p b (pink dashed line) for U a b ~ 0.1U a , U a = Ub, and fib = 0.75/x a and U a — 9 (b), 
U a = 11 (c), and U a = 13 (d). (We do not divide explicitly by zt because we set zt = 1). The phase diagram shows an SF 
phase and brown MI a lMIf,l, pink MI a 2MI&2 and red MI a 2MIi,l lobes; green MI a l, dark-green MI a 2 and a green-ochre MI 6 2 
regions; these are like their counterparts in Figs. [3(a) and (e). In addition we have a light-blue region MIf,l in which bosons of 
type b are in an MI phase with p b = 1 and bosons of type a are superfluid. In the second row we show the phase diagram (e), 
and plots versus p a of the order-parameters tj) a and ip b and the densities p a and pt for U a b = 0.3f7 a , Ub = U a , and pt = 0.75p a 
and U a = 9 (f), U a = 11 (g), and U a = 13 (h). The phase diagram shows the following: an SF phase; brown MI a l,MIf,l and 
red MI a 2MI;,l lobes; green MI a l, dark-green MI a 2 and light-blue MI(,1 regions; these are like their counterparts in (a). In 
addition we have an SF a phase in which bosons of type a are in an SF phase and the bosons density of type b are vanished. In 
the third row we show the phase diagram (i), and plots versus p a of the order-parameters xj) a and ipb and the densities p a and 
Pb for Uab = 0-7U a , Ub — U a , and p b = 0.75p a and U a = 9 (j), U a = 11 (k), and U a = 13 (1). The phase diagram shows the 
following: an SF a phase; MI a l and MI a 2 regions; these are like their counterparts in (a) and (e) in which the bosons density 
of type b are vanished. 



be possible to obtain them in time-of-flight measure- 
mentsi. Three-dimensional transforms of the shell struc- 
ture can be obtained, but they are not easy to visualize; 
therefore, we present the one-dimensional Fourier trans- 
forms of p a (X,Y = 0,Z = 0) and p b (X,Y = 0,Z = 0) 
with respect to X. The moduli of these transforms, 
namely, |p a ,fe x l (green lines), and |pb,fc x | (pink lines) of 
the profiles in Figs. [Til] (a)- (f) are plotted, respectively, in 
Figs. 1111 (a)-(f) versus the wave vector fcx/i"- The prin- 
cipal peaks in these transforms occur at kx = (or 2ir); 
these are associated with the spatially uniform MI and 
SF phases . In an infinite system with no confining po- 
tential, these are the only peaks; however, the quadratic 
confining potential leads to shells of MI and SF phases 
(see below); this shell structure leads to the subsidiary 
peaks that appear in Figs.[UJ(a) - (f) away from kx = 0, 
and 2tt. 

We can also obtain order-parameter-profile plots; these 



are shown in Figs.fTSTa) - (f) for the same parameter val- 
ues as in Figs. [TUT a) -(f), respectively; in these plots ip a is 
indicated by a red dashed line and ipb by a blue line. We 
can also obtain the one-dimensional Fourier transforms 
of ip a (X, Y = 0, Z = 0) and ip h {X, Y = 0, Z = 0) with 
respect to X. The moduli of these transforms, namely, 
iVVfcx I ( re d lines) and \ipb,kx I (blue lines) of the profiles 
Figs. [PHa) - (f) are plotted, respectively, in Figs. [TUTa) 
- (f ) versus the wave vector kx /tt . Again, the princi- 
pal peaks in these transforms occur at kx — (or 2tt); 
but subsidiary peaks occur because of the shell structure 
imposed by the confining potential. 

In Fig. [IJJa) - (d) we show representative density and 
order-parameter profiles, for p a and ip a , and the mod- 
uli of their Fourier transforms for U a b = 2.22U a , Ub — 
0.65U a , U a = 13, fib = 0-S^a,V T = 0.008, and L = 64. 
With this set of parameter values pb and ipb vanish. 

The density and order-parameter profiles of Figs. [TU] 
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Ub = 0.9U a , and /if, = 0.9/i a and U a = 9 (f), U a = 11 (g), and £/ a = 13 (h). The phase diagram shows the following: an SF 
phase and SF a ; brown MI a l,MIj,l, red MI a 2MIf,l and pink MI a 2MI 6 2 lobes; green MI a l, dark-green MI a 2 and light-blue MI 6 1 
regions; these are like their counterparts in Figs. [8] (a) and (e). 



and Figs. [T2l lead to the shell structure that we describe 
below. For specificity, consider Fig. fTOT d) and Fig. lT^T d'). 
These plots show that along the line Y = Z = 0, in the 
region from X = to \X\ = 10 the system has an MI 
phase for both types of bosons with p b = p a = 2 and 
ipa = ipb = 0; in the regions 10 < X < 28 and -28 < 
X < —10), the system displays an SF phase for both 
types of bosons with pb > p a and slightly ip a > ip b > 0; 
in the intervals 28 < X < 50 and -50 < X < -28 
bosons of type a are in an MI phase with p a = 1; when 
28 < X < 34 or -34 < X < -28 bosons of type b are 
in an MI phase with p b = 2; in the regions 34 < X < 42 
and —42 < X < —34 bosons of type b are in an SF 
phase whereas those of type a are still in the MI phase 
with p a = 1; if 42 < X < 50 or -50 < X < -42, 
then the displays MI phase for both types of bosons with 
Pb = Pa = 1; at slightly larger values of \X\ the system 
moves into an SF phase with p a — Pb > and ipa = 
tpb > 0; at even larger values of \X\ the system moves 
into a very narrow MI phase for the both types of bosons 
with p a = pb — 0.5, such that the total density for the 
system is p = p a + pb = 1 and %jj a — ipi — 0; as \X\ 
increases further, the system displays a very narrow SF 
region until it enters a small region in which the boson 
density vanishes for the both types of bosons. 

From the profiles in Figs. [TU] and Figs. [T2l fa) - (f) it is 
clear that the precise sequence and types of MI and SF 



shells depends on the parameters in the BH model ([3]) for 
two species of bosons. These SF and MI shells appear 
as annuli in a two-dimensional (2D) planar section Vz 
through the 3D lattice, at a vertical distance Z from the 
center as shown, for Vz=0: m Figs. [13 (a) and (b), for 
bosons of types a and b, respectively for the profiles of 
Figs.fTUTd) and [T2T dV Similar annular structures for the 
profiles of Figs. [10] (e) and [12] (e) are given in Figs. [15] 
(c) and (d) . We do not display the annular structures for 
the other profiles in Figs. [TU] and Figs. [T^] 

For any 2D planar section Vz we can calculate inte- 
grated, in-trap density profiles such as N m (Z), the num- 
ber of bosons in the p = m MI annuli, as we discussed 
for the BH model with one species of bosons. Here m 
is an integer; we concentrate on m = 1 or m = 2. We 
can also calculate the remaining number of bosons, e.g., 
N r am (Z) = N aT - N am (Z) or N^Z) = N bT - N bm (Z). 
For the parameter values of Figs. [TS] (a) and (b), illustra- 
tive integrated, in-trap density profiles are plotted ver- 
sus Z in Figs. [15] (e) and (f), respectively. These in-trap 
profiles show the total number of bosons for type a and 
b, N a T and N b T (blue full lines), the number of bosons in 
MI2 and Mil regions, N a2 , N b2 (white line in (e) and (f)) 
and N ai , N bl (green dash line in (e) and (f)), respectively, 
the numbers of bosons in MI( a ,b) 1/2 (red dash line in 
(e) and (f)) regions, [N aT — N a2 , N bT — N b2 ] (black full 
line in (e) and (f)), and [N aT — iV ai ,N bT — N bl ] (pink 
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FIG. 10: (Color online) Plots of p a % (green dashed line) and p bi (pink line) versus X, along the Y — Z = line, with 
Ma = Pb = 30, V T /(zt) = 0.008, and (a) U a /(zt) = 8 , U b = 0.9U a , and U ab = 0.6U a (b) U a /(zt) = 10 , U h = 0.9U a , 
and U ab = 0.6U a , (c) U a /(zt) = 11 , U„ = 0.9U a , and U ab = 0.6U a , (d)Z7«/(zt) = 12 , U b = 0.9U a , and U ab = 0.6U a , (e) 
J7 a /(^) = 13 , U b = 0.9U a , and U ab = 0.6U a , and (f) (7 n /(^) = 14 , {7 fe = 0.9Z7 a , and U ab = 0.6f7 a . 
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FIG. 11: (Color online) The moduli of the one-dimensional Fourier transforms, namely, \p( a . b )k\, of the plots of P( a , b )i in 
Figs. fTOT a 1 )- (f) are plotted, respectively, in (a) - (f) here versus the wave vector kx/~K. 
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FIG. 12: (Color online) Plots of Plots of ip a i (green dashed line) and ipw (pink line) versus X, along the Y — Z = line, 
with fi a = n b = 30, Vr/(zt) = 0.008, and (a) U a /(zt) = 8 , U b = 0.9U a , and U ab = 0.6U a (b) U a /{zt) = 10 , U b = 0.917a , 
and U ab = 0.6U a , (c) U a /(zt) = 11 , U b = O.W a , and U ab = O.W a , (d)U a /(zt) = 12 , U b = 0.9U a , and U ab = 0.6Z7«, (e) 
U a /(zt) = 13 , U b = 0.9!7 a , and j7 a6 = 0.6f/ a , and (f) U tt /(si) = 14 , U b = 0.9U a , and U ab = 0.6t/ a . 




FIG. 13: (Color online) The moduli of the one-dimensional Fourier transforms, namely, \ip( a , b )k\, of the plots of ip( a ,b)i in 
Figs. I12f a~) - (f) are plotted, respectively, in (a) - (f) here versus the wave vector hx/tr. 
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FIG. 14: (Color online) Plots of (a) p a ,i (b) V<M versus X, along the Y = Z = line, with U ab = 2.22U a , U b = 0.65C/ a , f7„ = 
13, /Xi, = 0.8/j.a, Vt = 0.008, and L — 64 (with this set of parameter values pt and ipb vanish); corresponding plots of the moduli 
of the one-dimensional Fourier transforms, namely, (c) |p a ,fc|, and (d) \i/) a ,k\ versus the wave vector kx- 




FIG. 15: (Color online) SF (black), MI2 (white), Mil (green), MI(a,b) 1/2 (red) annuli in the 2D planar section Vz=o (see text) 
with U a /(zt) = 13, Ut = 0.9U a , Uab = 0.6U a , Vr/(zt) — 0.008 and p a = jib = 30 (a) type a (b) type b; for these parameter 
values, the integrated, in-trap density profiles are plotted versus Z in (e) and (f), respectively. These in-trap profiles show the 
total number of bosons for type a and 6; N aT ,Nb T (blue full lines), the number of bosons in MI2 and Mil regions, N a2 , Nb 2 
(white line in (e) and (f)) and 7V ai , Nb 1 (green dash line in (e) and (f)), respectively, the numbers of bosons in MI(a,b) 1/2 
(red dash line in (e) and (f)) regions, [N aT ~ N a2 , Nb T — Nb 2 ] (black full line in (e) and (f)), and [N aT — 7V ai ,Nt T — N bl ] 
(pink annulus in (e) and (f)). The annular structures for the profiles of Figs. [101 (e) and [12] (e) are given in (c) and (d) and 
corresponding integrated, in-trap density profiles in (g) and (h), respectively. The outermost gray regions in (a) - (d) contain 
no bosons. 
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annulus) in (e) and (f). The outermost gray regions con- 
tain no bosons. Integrated, in-trap density profiles for 
the planar sections in Figs. [15] (c) and (d) are shown in 
Figs. H51 (g) and (h), respectively. 



C. Results for the Spin-1 Bose-Hubbard Model 

With the order parameters that we have defined in 
Eq. [10] we can, first, obtain phase diagrams for the spin- 
1 BH model for various values of Uq and ETjj; we re- 
fer the reader to our earlier study 18 for such phase di- 
agrams that include polar and ferromagnetic SF phases. 
Here we use the inhomogeneous MF theory, which we 
have developed above for the spin-1 BH model, to ob- 
tain some illustrative results for order-parameter profiles 
in a representative case that has a polar superfluid. In 
particular, we consider a simple-cubic lattice with 70 3 
sites, p/E r = 1, V T /E r = 0.001, V /E r = 14.5, and 
U2/U0 = 0.03, the 2D planar section V z for z = is 
plotted in Fig. (|16l) fa). The radial variations of the total 
on-site density of bosons pi — J2a Pi,a an d total on-site 
superfluid density p\ = J2 a Pt.a = J2a I ^,<r | 2 , as well as 
the individual component of superfluid density | | 2 
are given in Figs. <| 16|) (b>) and (c), respectively. From 
Fig. [ToTb) it is evident that this system has two well de- 
veloped MI (p — 1 and 2) shells, which are represented 
as regions with black and red respectively. The most im- 
portant result of model ([3]^ is that the superfluid phase 
is polar for U2 > and, according to the symmetry con- 
sideration, within our MF theory, the superfluid order 
parameters take one of the two possible set of values; 
ip±i ^ 0, Vo = or V±i = 0, + 0. Figure CEl])(c) 
yields tp±l > and i/'o = in the superfluid phase con- 
firming the polar nature of the phase. Another important 
feature is that pi ±i 7^ Pi ; o m the polar superfluid phase. 
This leads to N±i(z) 7^ Nq(z) where N a (z) is the total 
number of bosons with a spin a in the 2D planar sec- 
tion V z and is plotted in Fig. [ToTd) versus z. Thus the 
determination of N a (z) experimentally can reveal these 
features and thus can be used to confirm the polar nature 
of the superfluid phase in spin-1 bosons in optical lattice. 

In Fig. [T7] we show moduli of the one-dimensional 
Fourier transforms of the density and order-parameter 
profiles in Fig. Ql)] (c). It would be interesting to see if 
such patterns can be obtained via time-of-flight measure- 
ments. 

In Fig. [TH] we show a representative plot of the analog 
of Fig. [5] for the spin-1 BH model with parameter val- 
ues as in Fig. QIdJ thus, there are two well-developed MI 
shells. Here N^{p) denotes the total number of bosons, 
at density p, and with a = ±1; similarly, N®(p) is the to- 
tal number of bosons, at density p, and with a = 0; and is 
Ni,(p) is the total number of bosons at density p. For the 
peak in ./Vf,(p), in the vicinity of p = 1, only bosons with 
a = ±1 contribute; but, for the one near p = 2, all three 
components contribute equally. This result, which is also 
implicit Fig. [in] (c), should be verifiable in atomic-clock- 



shift experiments of the type that have been carried out 
for spinless bosons^. 

We have noted in earlier worki£ that our mean-field 
theory does not account for order parameters that distin- 
guish between different spin orderings, which have been 
studied^ in the limit Uq — > 00, in the MI phases in spin- 
1 BH models. The exploration of such orderings lies be- 
yond the scope of the present study. 
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FIG. 16: (Color online) (a) SF (white) and MI regions 
[p — 2(red) and p = 1 (black)] annuli formed in the 2D section 
though the 3D optical lattice at the center and (b) the corre- 
sponding radial variation of density pi and superfluid density 
pi. (c) radial variation of pi a and p\.& and (d) N a {z) versus 
z for cr = ±1,0. Note that N-i(z) =Ni(z). N(z) = £ CT N a . 
Here the optical lattice parameters are taken to be p/E r — 1, 
Vr/Er = 0.001, Vo/Er = 14.5 and U 2 /U = 0.03. 



IV. CONCLUSIONS 

We have carried out a comprehensive study of Mott in- 
sulator and superfluid shells in Bose-Hubbard models for 
bosons in optical lattices with harmonic traps by using an 
intuitively appealing inhomogeneous mean field theory 
that has been used earlier to understand the Bose-glass 
phased. Our inhomogeneous mean-field theory quanti- 
tatively agrees with QMC simulations. Furthermore, it is 
numerically less intensive than QMC simulations; thus, 
we are able to perform calculation on experimentally re- 
alistic, large 3D systems and explore a wide range of pa- 
rameter values. We can calculate in-trap density profiles 
that agree qualitatively with experiments^; and we show 
how to obtain the phase diagram of the homogeneous 
Bose Hubbard model from such in-trap density profiles. 
Our results are also of direct relevance to recent atomic- 
clock-shift experiments^ as we have described above. Fi- 



16 




15 



(d) 



10 



0.5 1 1.5 2 

k /it 



.^Jl 

0.5 1 1.5 2 

k In 



FIG. 17: (Color online) Illustrative plots of the moduli of the one-dimensional Fourier transforms (a) IV'o.fcli (b) IV'i.kli ( c ) 
|p ,fe|, and (d) |pi,fe| for the profiles in Fig. [16] (c). 
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FIG. 18: (Color online) Representative plots versus the den- 
sity p of N^(p) (red circles), the total number of bosons, at 
density p and with a = ±1, N®(p) (black squares), the total 
number of bosons, at density p and with a = 0, and their 
sum N b (p) = Nf(p) + N°(p) + N^(p) (blue triangles); here 
all parameters are as in Fig. 1161 



nally we have generalized our inhomogeneous mean-field 
theory to BH models with two species of bosons or a spin- 
1 BH model with harmonic traps. With two species of 
bosons we obtain rich phase diagrams with a variety of 
SF and MI phases and associated shells, when we include 
a quadratic confining potential; we also obtain in-trap 
density distributions that show plateaux as in the single- 
species case. For the spin-1 BH model we show, in a rep- 
resentative case, that the system can display alternating 
shells of polar SF— and MI phases; and we make inter- 
esting predictions for atomic-clock-shift experiments. We 
hope our results will stimulate more experiments on such 
systems of bosons in optical lattices. Our inhomogeneous 
mean-field theory can also be generalized to study the ex- 
tended Bose-Hubbard model as we report elsewhere^. 

Though other groups^— have studied such shell 
structure theoretically, they have not obtained the quan- 
titative agreement with quantum Monte Carlo (QMC) 
simulations 2 ^ that we obtain, except in one dimension 3 ^. 
Furthermore, there have been some investigations of the 
BH model with a harmonic trap potential; these use 
mean-field theor y 13 ' 28 and, in addition, a local-density 
approximation (LDA) , which assumes that the properties 
of a system with finite confining potential at a particular 
location are identical to those of a uniform system with 
the value of the local chemical potential at that location. 
This approximation leads to a decoupling of each site 
from its neighbor; it is equivalent to assuming 0; = ipi in 
Eq. O and a minimization of the ground-state energy for 
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each site separately. In our inhomogeneous mean-field 
theory, we do not make this additional LDA assumption; 
and the minimization of the ground state energy is done 
over the entire set of ipi. If we compare these two ap- 
proaches for the single-species BH model, we find that 
the difference is negligible in SF regions, but discrepan- 
cies exist at SF-MI interfaces; this has been reported in 
other models^ also. 
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